Condensate statistics in interacting Bose gases: exact results 
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Recently, a Quantum Monte Carlo method alternative to the Path Integral Monte Carlo method 
was developed for the numerical solution of the iV-boson problem; it is based on the stochastic 
evolution of classical fields. Here we apply it to obtain exact results for the occupation statistics of 
the condensate mode in a weakly interacting trapped one-dimensional Bose gas. The temperature 
is varied across the critical region down to temperatures lower than the trap level spacing. We 
verify that the number-conserving Bogoliubov theory gives accurate predictions provided that the 
non-condensed fraction is small. 
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The first achievement of Bose-Einstein condensation in 
weakly interacting atomic gases in 1995 has renewed the 
interest on basic aspects of the Bose-Einstein condensa- 
tion phase transition jij. In particular, an intense the- 
oretical activity has been recently devoted to the study 
of the occupation statistics of the condensate mode: al- 
though no experimental result is available yet, we expect 
that it would provide a crucial test for many-body theo- 
ries since, contrarily to more common one-body observ- 
ables like the density, it involves arbitrarily high order 
correlation function of the quantum Bose field. While 
there are now well established results for the ideal Bose 
gas (|], calculations for the weakly interacting Bose gas 
have been performed within the framework of mean-field 
approximation only |j, Q ||. The intermediate regime 
around the critical temperature where the mean-field the- 
ory fails has therefore been left unexplored. Some con- 
troversy is still open about the validity of the Bogoliubov 
approach even at temperatures much lower than the crit- 
ical temperature |6). In the absence of experimental re- 
sults, it is then very interesting to have exact theoretical 
results on the statistics of condensate occupation. 

There exists an exact analytical solution to the bosonic 
N body problem Q, but it is restricted to the spatially 
homogeneous one-dimensional case. From the side of nu- 
merics, the Quantum Monte Carlo method based on the 
Path Integral Monte Carlo technique is in principle able 
to give exact predictions for any observable of the gas || 
and was successfully used to calculate the mean conden- 
sate occupation || and the critical temperature [^0| . In 
the Path Integral formulation, however, the position rep- 
resentation is priviledged which makes the calculation of 
highly non-local observables like the condensate occupa- 
tion probabilities rather involved. 

In this paper, we use instead a recently developed 
Quantum Monte Carlo method based on the stochastic 
evolution of classical fields |ll], |l2) . This new method has 
a much broader range of applicability than the standard 
Path Integral Monte Carlo technique: it can be applied to 
bosonic systems with complex wavefunctions, and even 



to interacting Fermi systems |L3[ . As compared to the 
Positive- P method of quantum optics jl4j , it has the de- 
cisive advantage of having been proven to be convergent. 
In this paper, we perform the first non-trivial applica- 
tion of this new method, calculating for the first time the 
exact distribution function of the number of condensate 
particles in the presence of interactions, for temperatures 
across the Bose-Einstein condensation temperature. 

An ultracold trapped interacting Bose gas in D dimen- 
sions can be modelled in a second-quantization formalism 
by the Hamiltonian 

H = E |^"k«k + I £ dV *t (r )*t ( r )*( r )*( r ) 

k r 

+ Y / d VU cxt (r)¥(r)^(r); (1) 

r 

the spatial coordinate r runs on a discrete orthogonal lat- 
tice of J\f points with periodic boundary conditions; V is 
the total volume of the quantization box and dV = V/JV 
is the volume of the unit cell of the lattice. C/ GX t(r) is the 
external trapping potential, m is the atomic mass and 
interactions are modeled by a two-body discrete delta 
potential with a coupling constant go. The field opera- 
tors ^(r) satisfy the usual Bose commutation relations 
[*(r), ^(r')] = S r ^ r i/dV and can be expanded on plane 
waves according to ^(r) = ^2 k dke lkr / W with k re- 
stricted to the first Brillouin zone of the reciprocal lattice. 
In order for the discrete model to correctly reproduce the 
underlying continuous field theory, the grid spacing must 
be smaller than the macroscopic length scales of the sys- 
tem like the thermal wavelength and the healing length. 

In this paper, we assume that the gas is at thermal 
equilibrium at temperature T in the canonical ensemble 
so that the unnormalized density operator is p eq (P) = 
e~P n with (3 = 1/ksT. It is a well-known fact of quan- 
tum statistical mechanics that this density operator can 
be obtained as the result of an imaginary-time evolution: 
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during a "time" interval r = — > /3 starting from the 
infinite temperature value p e q{j = 0) = In- We have 
recently shown |f2y that the solution of the imaginary- 
time evolution (0) can be written exactly as a statistical 
average of Hartree operators of the form 



g = |JV:&) (N: 2 1 



(3) 



where, in both the bra and the ket, all N atoms share 
the same (not necessarily normalized) wave functions <p a 
(a = 1,2). For the model Hamiltonian ([j]) here consid- 
ered, this holds if each cf> a evolves according to the Ito 
stochastic differential equations: 



dr 

.M = -y 



_ g (N-l) Y^dV |0 a (r') 
2 ll^ll 4 
with a noise term given by 



a (r) + dB Q (r) (4) 



dB a {r) 



2V 



Qa [4>a(l 



y^(k.r+J a (kj) 
k>0 



2.C.)] 



(5) 

where the projector Q a projects orthogonally to <j) a , the 
index k is restricted to a half space and the 9 a (k) are in- 
dependent random angles uniformly distributed in [0, 2n]. 

Starting from this very simple but exact stochastic for- 
mulation, a Monte Carlo code was written in order to 
numerically solve the stochastic differential equation 
As a first step, a sampling of the infinite temperature 
density operator has to be performed in terms of a fi- 
nite number of random wave functions <j>W. Since the 
effective contributions of the different realizations to the 
final averages can be enormously different, the statistics 
of the Monte Carlo results can be improved by using an 
importance sampling technique M so to avoid wasting 
computational time. This is done using the identities 



1 N = V(t>\N :c}>){N :<t>\ 



p [<p]v<p 



\N:<f>)(N: 



(6) 

where the integration is performed over the unit sphere 
\\(j>\\ 2 — {(j)\4>) = J2 T dV \(j>\ 2 = 1 and where the a priori 
distribution function Po[0] can be freely chosen in order 
to maximize the efficiency of the calculation. For the 
numerical calculations here reported, the following Pq[4>] 
has been used 



Pi 



2 A' 



(7) 



since it joined the possibility of a simple sampling with 
a fast convergence. In this expression, hap is the 

2 

Gross-Pitaevskii Hamiltonian hop = #— + J7 cxt (r) + 
N go |(/)q (r)| 2 — /i, 4>q is the wave function which minimises 



the Gross-Pitaevskii energy functional and \i is the cor- 
responding chemical potential. With respect to the ideal 
Bose gas distribution function previously used [Q, the 
present form (Q) for Po[(f>] has the advantage of taking 
into account the fact that the condensate mode can be 
strongly modified by interactions. More details on the 
actual sampling of Pq [</>] can be found in |12| . 



Each realization 



/.(*) 



is then let evolve according to 



the stochastic evolution in imaginary-time (|4|) from its 
t = value <f>® 2 = <f>® to the inverse temperature of 
interest t — [3. The expectation values of any observable 
at temperature T can then be calculated as averages over 
the Monte Carlo realizations; e.g., the partition function 
Tr [p] is given by Tr[p] = (4>2\4>i) N ■ I n particular, the 
condensate wavefunction 0BEc( r ) is the eigenvector of 
the one-body density matrix 

1 



*t(r')#(r) 



Tr[p] 



\N-l 



(8) 



corresponding to the largest eigenvalue, that is to the 
largest mean occupation number Jl5|] . The complete 
probability distribution P(N ) for the occupation statis- 
tics of the condensate mode is obtained via the expression 
P(No) = Tr[pN p] Tr[p] _1 where Vn projects onto the 
subspace in which the condensate mode contains exactly 
No atoms; in terms of the ansatz (|J), the expectation 
value of the projector can be written as 

Tr[V NoP ] = ^ No)] (c| Cl )"o <#|#)"-*o ( 9 ) 

where we have split (f> a (r) as c Q 0BEc( r ) + 0a( r ) with 
4>^ (r) orthogonal to the condensate wavefunction. 

In fig.Q-||we have summarized the results of the Monte 
Carlo calculations for a one-dimensional, harmonically 
trapped gas with repulsive interactions. In fig.|l|, we 
have plotted the condensate wavefunction density pro- 
file |</>bec(£)| 2 f° r different values of the temperature, 
see solid lines. As the temperature decreases, the en- 
hanced atomic density at the center of the cloud gives 
a stronger repulsion among the condensate atoms and 
therefore a wider condensate wavefunction. In this low 
temperature regime, we have found a good agreement be- 
tween the Monte Carlo results and the Bogoliubov predic- 
tion for the condensate wavefunction including the first 
correction to the Gross-Pitaevskii equation due to the 
presence of non-condensed atoms Q, see dashed lines 
in fig.§ This was expected, since the numerical exam- 
ples in this paper are in the weakly interacting regime 
n£ w 15 3> 1, where n is the density at the trap cen- 
ter and £ = (fi 2 /2mgon) 1 / 2 is the corresponding heal- 
ing length. For high temperatures, the wavefunction 
0BEc(aO tends to the harmonic oscillator ground state. 

The wavefunction 4>bec(x) can then be used to deter- 
mine the occupation statistics of the condensate mode 
via (|^). The signature of Bose condensation is appar- 
ent in fig.pl: the occupation statistics of the condensate 
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mode at high temperatures has the same shape as a 
non-interacting, thermally occupied mode, its maximum 
value being at N — 0. For decreasing T, P(N ) rad- 
ically changes its shape; at low T, its maximum is at 
a non-vanishing value of Nq and its shape is strongly 
asymmetric with a longer tail going towards the lower 
Nq values, as already noticed in the mean-held approach 
of p|. This transition ressembles the one occurring in a 
laser cavity for a pumping rate which goes from below to 
above threshold p7[ . 

For even lower values of the temperature, see fig||, 
most of the atoms are in the condensate and the probabil- 
ity distribution P(Nq) tends to concentrate around val- 
ues of No close to N. However, interactions prevent the 
atomic sample from being totally Bose condensed even 
at extremely low temperatures. Moreover, the probabil- 
ity of having an odd number of non-condensed atoms is 
strongly reduced with respect to the probability of hav- 
ing an even number (see the strong oscillations in the 
ksT = 0.4 Hlo curve of hg.^). A quantitative interpreta- 
tion of this effect in terms of the Bogoliubov approxima- 
tion will be given in the following part of the paper. 

We now compare the exact results with the prediction 
of the number- conserving Bogoliubov approximation Jji], 
|l8f . The Bogoliubov Hamitonian has the form 



with 



£ 



H Bos = ^dV (At, A 



hap + QoNg\<f>o\ Qo 



A 

At 



(10) 



QoNg<ftQ& 



-Q* Ngr 2 Q ~h GP - Q* Ng |</> | 2 Q* 



The matrix 77 is defined according to 

fl 
^= 0-1 



(11) 

(12) 



and the projector Qo projects orthogonally to the Gross- 
Pitaevskii ground state 4>o- The operators A are defined 
according to (A) r = Aj^^r) where $x( r ) is the projec- 
tion of the field operator orthogonally to (j> and where 
A 1 N : 0o ) = I N - 1 : (f> ) if N > and otherwise. Ne- 
glecting the possibility of completely emptying the con- 
densate mode, Ao and Aj can be shown to commute. Un- 
der this approximation, the A are bosonic operators. The 
probability H(5N) of having 6N non-condensed atoms 
can be calculated within the Bogoliubov approximation 
by means of the characteristic function H 



/(<?) =J2 n i 6N ) elSNe - Tr 



6N 



18SN ± -/3Hb os 

z 



(13) 



when the 



with SN = dV At • A and Z = Tr [e 
probability of emptying the condensate mode is not to- 
tally negligible, the Bogoliubov prediction for IL(SN) has 



to be truncated to the physically relevant 5N < N values 
and then renormalized to 1 . The details of the calcula- 
tion will be given elsewhere, we go here to the result 



2(JV-1) 

n 

3=1 



2X< 



l + Aj-e^l-Aj) 



1/2 



(14) 



in which the Xj are the eigenvalues of M. = rytanh(/3£/2); 
from such an expression, the occupation probabilities 
IL(8N) are immediately determined by means of an in- 
verse Fourier transform. The results are shown as dashed 
lines in figs.^ and ^ where they are compared to the ex- 
act results from Monte Carlo simulations; the agreement 
is excellent provided the non-condensed fraction is small, 
i.e. at sufficiently low temperatures. 
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FIG. 1: Solid lines: exact results for the condensate wave- 
function density profile |</>BEc(a:)| 2 for N = 125 atoms in 
a one- dimensional harmonic trap at different temperatures 
k B T/huj = 1 (a), 10 (b), 20 (c), 50 (d). lo is the trap frequency, 
|</>bec| 2 is in units of the inverse of <Zho = \J h/mui. The cou- 
pling constant is go = 0.08 huj aha- Comparison with Bogoli- 
ubov predictions (dashed lines) in (a,b,c) and with harmonic 
oscillator ground state (diamonds) in (d) . In (a) the two lines 
are indistinguishable on the scale of the figure. Monte Carlo 
calculations were performed using 3 ■ 10 4 realisations on a 
j\f =16 (a), 64 (b,c) or 128 point (d) grid. 



In order to understand the oscillating behaviour shown 
by the lowest temperature curve of fig.[3|, we calculate the 
ratio of the probability n o dd of having an odd value of 
5N to the probability n cV on of having an even value: 



nodd = /(Q) - /(tt) i-r 

Ileven /(0) + /(tt) 1 + R 



(15) 



where R = ]}. X 1 / 2 = det[M}^ 2 = U m tanh(/3e m /2) and 
m runs over the M— 1 eigenenergies e rn of the Bogoliubov 
spectrum. In order to have oscillations in H(5N), the ra- 
tio (|l^) must be small as compared to 1 which requires a 
temperature smaller than the energy of the lowest Bogoli- 
ubov mode. These oscillations are therefore a property 
of the ground state of the system. In the Bogoliubov 
approximation, the Hamiltonian is quadratic in the field 
operators so that its ground state is a squeezed vacuum, 
which indeed contains only even values of SN, a well 
known fact of quantum optics UM- From a condensed 
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FIG. 2: Solid lines: exact results for the distribution function 
P(No) of the number of condensate atoms for decreasing tem- 
peratures from left to right: k B T/huj — 50 (a), 33 (b), 20 (c), 
10 (d), 5 (e). Same parameters as in fig|l[ For curves (c,d,e), 
comparison with Bogoliubov predictions (dashed lines); for 
curve (a), comparison with the ideal gas (diamonds). The 
Monte Carlo calculations were performed using 3 • 10 4 reali- 
sations on a grid of TV = 128 (a,b) or N = 64 points (c,d,e). 



value of Nq. At temperatures below the trap oscillation 
frequency, configurations with an odd number of non- 
condensed particles are strongly suppressed, which we 
interpret successfully within the Bogoliubov approxima- 
tion. Possible extensions of this work are (i) the deter- 
mination of critical temperature shift for an interacting 
Bose gas in the limit of a vanishing scattering length, still 
subject of some controversy p2| , (ii) exact calculations 
for finite temperature Bose gases with vortices in rotat- 
ing traps, and (iii) the determination of the BCS critical 
temperature in two-component Fermi gases. 
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FIG. 3: Solid lines: exact results (10 5 realizations on a 
N — 16 point grid) for the distribution function H(SN) = 
P(No = N — SN) of the number of non-condensed atoms at 
low temperatures fesT/Tiw = 0.4 (a) and 1 (b). Same param- 
eters as in fig.[j] and fig.^. Dashed lines: Bogoliubov predic- 
tions. For kaT/fko = 0.4, the two lines are indistinguishable 
on the scale of the figure. 



matter point of view, several existing variational ansatz 
for the ground state corresponding to a condensation of 
pairs also exhibit this property [^, [H]. Physically, this 
is ultimately related to the fact that the leading interac- 
tion process changing the condensate occupation at zero 
temperature, that is the scattering of two condensate par- 
ticles into two non-condensed modes and vice-versa, does 
not change the parity of SN, as it is apparent in the Bo- 
goliubov Hamiltonian (]Tc|) . 

In conclusion, we have determined with a newly de- 
veloped stochastic field Quantum Monte Carlo method 
the exact distribution function of the number Nq of con- 
densate particles in a one-dimensional weakly interact- 
ing trapped gas. The signature of Bose condensation is 
the appearance of a finite value for the most probable 
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